
function J=J_bvg_thetas(U1,U0,R,theta,sgi_nodes,DS,DB)

N=length(DS)/2;

J=(R*theta(:,2)).*U1.*(DS(N+(1:N),1)*sgi_nodes'+DB(N+(1:N),1))-...
    (R*theta(:,1)).*U0.*(DS(1:N,1)*sgi_nodes'+DB(1:N,1));
